{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BenchmarkTools, Test\n",
    "using Images\n",
    "using OffsetArrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Julia 的下标索引规则\n",
    "\n",
    "在对多维数组进行下标索引 `A[I_1, I_2, ..., I_n]` 时，每一个 `I_k` 都可以是下面三种情况的任意一种：\n",
    "\n",
    "* 标量，例如 `A[1]` `A[1, 2]`\n",
    "* 一组标量组成的矩阵（或可以像矩阵一样工作的结构），例如 `A[[1, 3]]`, `A[1:3]`\n",
    "* 可以转换成第二种情况的数据，例如 `A[:, :]`\n",
    "\n",
    "针对这几种索引的类型，Julia 有三条索引的规则：\n",
    "\n",
    "* 坐标索引 (Cartesian Indexing)，例如 `A[1, 2, 3]`\n",
    "* 线性索引 (Linear Indexing)，例如 `A[5]`\n",
    "* 特殊情况：是否可以省略额外的 `:`\n",
    "\n",
    "下面先介绍这三种规则下各种类型的索引是如何使用的"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×4×2 Array{Int64,3}:\n",
       "[:, :, 1] =\n",
       " 1  3  5  7\n",
       " 2  4  6  8\n",
       "\n",
       "[:, :, 2] =\n",
       "  9  11  13  15\n",
       " 10  12  14  16"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 简单起见，先构造一个比较小的数组方便观察\n",
    "A = reshape(collect(1:16), 2, 4, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. 坐标索引\n",
    "\n",
    "常见的 `A[1, 2]` 这种通过指定数组在每一个维度的位置来进行索引的方式，称之为坐标索引。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "仅用标量进行坐标索引，最后得到的也是标量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1, 2, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "若索引中某一项改为第二类索引，即用一组标量组成的矩阵，那么最后得到的也是矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 11\n",
       " 15"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1, [2, 4], 2] # 等价于 [A[1, 2, 2], A[1, 4, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Int64,2}:\n",
       " 11  15\n",
       " 12  16"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[[1, 2], [2, 4], 2] # 等价于 [A[1, 2, 2] A[1, 4, 2]; A[2, 2, 2] A[2, 4, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Int64,1}:\n",
       " 11\n",
       " 11\n",
       " 11"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 并不要求矩阵中每个元素都是唯一的\n",
    "A[1, [2, 2, 2], 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "类似地，矩阵也可以替换成各种 Range 类型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 11\n",
       " 15"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1, 2:2:4, 2] # collect(2:2:4) == [2, 4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在图像处理中，索引图像就是以这种方式工作的，例如"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGQAAABkBAMAAACCzIhnAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAD1BMVEUAAAD/AAAA/wAAAP////+B57DtAAAAAWJLR0QEj2jZUQAAAEtJREFUWMPtzjENwEAMBEFTCIVQeArhjynduXUVRdZs7Tm5Kl3pTic9qRAEWUB+9g6CIB+R6WFPIwiygUwPexpBkA1ketjTCIIsIC9G9Jgv/mYpQwAAAABJRU5ErkJggg==",
      "text/plain": [
       "5×5 Array{RGB{Float64},2} with eltype RGB{Float64}:\n",
       " RGB{Float64}(0.0,0.0,0.0)  …  RGB{Float64}(1.0,1.0,1.0)\n",
       " RGB{Float64}(1.0,0.0,0.0)     RGB{Float64}(0.0,0.0,0.0)\n",
       " RGB{Float64}(0.0,1.0,0.0)     RGB{Float64}(1.0,0.0,0.0)\n",
       " RGB{Float64}(0.0,0.0,1.0)     RGB{Float64}(0.0,1.0,0.0)\n",
       " RGB{Float64}(1.0,1.0,1.0)     RGB{Float64}(0.0,0.0,1.0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 调色盘 = [黑，红，绿，蓝，白]\n",
    "palette = [RGB(0.0, 0.0, 0.0), RGB(1.0, 0.0, 0.0), RGB(0.0, 1.0, 0.0), RGB(0.0, 0.0, 1.0), RGB(1.0, 1.0, 1.0)]\n",
    "\n",
    "# idx中的每一项表示它对应的调色盘中的索引\n",
    "idx = [\n",
    "    1 2 3 4 5; # 黑，红，绿，蓝，白\n",
    "    2 3 4 5 1; # 红，绿，蓝，白，黑\n",
    "    3 4 5 1 2;\n",
    "    4 5 1 2 3;\n",
    "    5 1 2 3 4\n",
    "]\n",
    "\n",
    "palette[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有两种比较典型的第三类索引：`:` (`Colon`) 与 二值矩阵\n",
    "\n",
    "* 用 `:` 作下标可以作为占位符使用，表示选取对应维度的全部下标\n",
    "* 用二值矩阵作为下标，则表示选取其中所有值为 `true` 的量\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 11\n",
       " 12"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[:, 2, 2] # 等价于 A[axes(A, 1), 2, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×4 Array{Int64,2}:\n",
       "  9  11  13  15\n",
       " 10  12  14  16"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[:, :, 2] # 等价于 A[axes(A, 1), axes(A, 2), 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注：不同于numpy，用 `:` 作为下标，只选取当前的维度而不会扩散到其他维度，因此在这个例子中，`A[1, :]` 或者 `A[:, 1]` 这种形式都是不可行的，需要将第三个维度的坐标加上才行"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×2 Array{Int64,2}:\n",
       " 1   9\n",
       " 3  11\n",
       " 5  13\n",
       " 7  15"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1, :, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "用逻辑值作为下标可以简单快速地把矩阵的内容筛选出来，例如："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×4×2 BitArray{3}:\n",
       "[:, :, 1] =\n",
       " 0  0  0  0\n",
       " 0  0  0  0\n",
       "\n",
       "[:, :, 2] =\n",
       " 0  1  1  1\n",
       " 0  1  1  1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = A .> 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Int64,1}:\n",
       " 11\n",
       " 12\n",
       " 13\n",
       " 14\n",
       " 15\n",
       " 16"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[I]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注：用逻辑值作为下标，则必须保证索引矩阵 `I` 的尺寸与 `A` 的尺寸一致"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. 线性索引"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算机在读取连续的内存块时性能要优于读取不连续的内存块，因此保证内存连续是处理多维数组时必须使用的技巧\n",
    "\n",
    "如果矩阵是按列存储，而按行读取的话，那么内存是不连续的，转置矩阵就是这样一个例子"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.685 μs (2 allocations: 32.08 KiB)\n",
      "  6.112 μs (2 allocations: 32.08 KiB)\n"
     ]
    }
   ],
   "source": [
    "X = rand(64, 64);\n",
    "\n",
    "@btime collect($(X));\n",
    "@btime collect($(X')); # 将矩阵转置之后，内存不再是连续的了，因此索引的结果更慢"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当对多维数组进行索引时只传递了一个下标，那么就会采用线性索引。\n",
    "\n",
    "简单来说，**线性索引**就是将数组按照它在内存中的存储顺序展开成一列向量，然后在向量上进行对应的索引。\n",
    "\n",
    "使用线性索引的好处是它能够保证数据一定是按照内存连续的方式进行读取"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 在 Matlab 中，也可以这样做\n",
    "# 在 numpy 中，这种操作是会报错的\n",
    "\n",
    "A[15] # 等价于 A[1, 4, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 类似地，\n",
    "A[length(A)] == A[end]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "my_euclidean2 (generic function with 1 method)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 以计算两个数组的欧式距离为例 d = sqrt(∑(Xᵢ - Yᵢ))\n",
    "# 这是两种基本等价的迭代方式\n",
    "\n",
    "function my_euclidean1(X, Y)\n",
    "    retval = 0.0\n",
    "    @assert axes(X) == axes(Y)\n",
    "    \n",
    "    for i in 1:length(X)\n",
    "        retval += (X[i] - Y[i])^2\n",
    "    end\n",
    "    return sqrt(retval)\n",
    "end\n",
    "\n",
    "function my_euclidean2(X, Y)\n",
    "    retval = 0.0\n",
    "    @assert axes(X) == axes(Y)\n",
    "    \n",
    "    for (x, y) in zip(X, Y)\n",
    "        retval += (x - y)^2\n",
    "    end\n",
    "    return sqrt(retval)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  3.662 μs (0 allocations: 0 bytes)\n",
      "  3.578 μs (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "X = rand(64, 64);\n",
    "Y = rand(64, 64);\n",
    "\n",
    "@btime my_euclidean1($X, $Y);\n",
    "@btime my_euclidean2($X, $Y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16-element Array{Int64,1}:\n",
       "  1\n",
       "  2\n",
       "  3\n",
       "  4\n",
       "  5\n",
       "  6\n",
       "  7\n",
       "  8\n",
       "  9\n",
       " 10\n",
       " 11\n",
       " 12\n",
       " 13\n",
       " 14\n",
       " 15\n",
       " 16"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# : 扩散到其他维度只在线性索引时才能工作\n",
    "\n",
    "A[:] # 等价于 vec(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.3. 数组索引的一些特殊规则\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 若接下来的维度都只含有一个元素，则可以省略下标"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = reshape(collect(1:24), 2, 3, 1, 4, 1, 1, 1)\n",
    "\n",
    "A[1, 1, 1, 1] # 可以省略第五个下标"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 若接下来的索引全为1， 则可以超出数组的维数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = reshape(collect(1:24), 2, 3, 4)\n",
    "\n",
    "A[2, 3, 4, 1, 1, 1] # works, too"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. CartesianIndex 与 多维数组"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "线性索引的好处是它可以保证内存的连续性，但缺点是它丢失了高维矩阵的位置信息，因此在需要涉及到对高维矩阵的位置进行交互的时候，往往还是会回到坐标索引下去工作。\n",
    "\n",
    "滤波、卷积这类运算就是非常典型的一个例子\n",
    "\n",
    "以均值滤波为例，在其他语言中一般会对不同维度提供一个特定的实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10×2 Array{Float64,2}:\n",
       " 5.0  4.0\n",
       " 3.0  4.0\n",
       " 4.0  4.0\n",
       " 5.0  4.0\n",
       " 3.0  3.66667\n",
       " 3.0  3.33333\n",
       " 4.0  2.66667\n",
       " 1.0  2.0\n",
       " 1.0  1.33333\n",
       " 2.0  1.5"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mean_filter_1d(X::AbstractVector; window_size=3)\n",
    "    out = similar(X)\n",
    "    \n",
    "    δ = window_size ÷ 2\n",
    "    for i in 1:length(X)\n",
    "        # min, max确保下标不溢出\n",
    "        l = max(1, i-δ)\n",
    "        r = min(length(X), i+δ)\n",
    "        \n",
    "        # 计算一个局部的均值\n",
    "        out[i] = mean(X[l:r])\n",
    "    end\n",
    "    return out\n",
    "end\n",
    "\n",
    "X = Float64.(rand(1:5, 10))\n",
    "hcat(X, mean_filter_1d(X; window_size=3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×10 Array{Float64,2}:\n",
       " 3.0  1.0  2.0  5.0  1.0  1.75     2.16667  2.5      2.33333  2.0\n",
       " 1.0  2.0  4.0  1.0  1.0  2.5      2.66667  3.0      2.66667  2.5\n",
       " 4.0  4.0  3.0  5.0  2.0  3.33333  3.44444  3.11111  2.44444  1.83333\n",
       " 5.0  4.0  4.0  1.0  1.0  3.5      3.33333  2.88889  2.44444  2.16667\n",
       " 3.0  1.0  2.0  2.0  2.0  3.25     3.16667  2.33333  2.0      1.5"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mean_filter_2d(X::AbstractMatrix; window_size=3)\n",
    "    out = similar(X)\n",
    "    M, N = size(X)\n",
    "    \n",
    "    δ = window_size ÷ 2\n",
    "    for i in 1:M\n",
    "        il = max(1, i-δ)\n",
    "        ir = min(M, i+δ)\n",
    "        \n",
    "        for j in 1:N\n",
    "            jl = max(1, j-δ)\n",
    "            jr = min(N, j+δ)\n",
    "        \n",
    "            out[i, j] = mean(X[il:ir, jl:jr])\n",
    "        end\n",
    "    end\n",
    "    return out\n",
    "end\n",
    "\n",
    "X = Float64.(rand(1:5, 5, 5))\n",
    "hcat(X, mean_filter_2d(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia 的 `CartesianIndex` 对于这种情形提供了一个简洁且高效的解决方案，在给出新的实现之前，先介绍 `CartesianIndex` 的基本使用"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`CartesianIndex` 只是对 `Tuple` 类型的一个包装"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 2, 3)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I_t = (1, 2, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3×4 Array{Int64,3}:\n",
       "[:, :, 1] =\n",
       " 1  3  5\n",
       " 2  4  6\n",
       "\n",
       "[:, :, 2] =\n",
       " 7   9  11\n",
       " 8  10  12\n",
       "\n",
       "[:, :, 3] =\n",
       " 13  15  17\n",
       " 14  16  18\n",
       "\n",
       "[:, :, 4] =\n",
       " 19  21  23\n",
       " 20  22  24"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = reshape(collect(1:24), 2, 3, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "ename": "ArgumentError",
     "evalue": "ArgumentError: invalid index: (1, 2, 3) of type Tuple{Int64,Int64,Int64}",
     "output_type": "error",
     "traceback": [
      "ArgumentError: invalid index: (1, 2, 3) of type Tuple{Int64,Int64,Int64}",
      "",
      "Stacktrace:",
      " [1] to_index(::Tuple{Int64,Int64,Int64}) at ./indices.jl:297",
      " [2] to_index(::Array{Int64,3}, ::Tuple{Int64,Int64,Int64}) at ./indices.jl:274",
      " [3] to_indices at ./indices.jl:325 [inlined]",
      " [4] to_indices at ./indices.jl:322 [inlined]",
      " [5] getindex(::Array{Int64,3}, ::Tuple{Int64,Int64,Int64}) at ./abstractarray.jl:1060",
      " [6] top-level scope at In[27]:1"
     ]
    }
   ],
   "source": [
    "A[I_t] # Tuple 类型不能作为索引使用"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CartesianIndex(1, 2, 3)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = CartesianIndex(I_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[I] # 等价于 A[1, 2, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在对 `Tuple` 进行包装的基础上，`CartesianIndex` 上有一些额外的下标操作"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "比如说，用 `:` 来构造一组 `CartesianIndex` 的下标"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×2×3 CartesianIndices{3,Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}:\n",
       "[:, :, 1] =\n",
       " CartesianIndex(1, 2, 2)  CartesianIndex(1, 3, 2)\n",
       "\n",
       "[:, :, 2] =\n",
       " CartesianIndex(1, 2, 3)  CartesianIndex(1, 3, 3)\n",
       "\n",
       "[:, :, 3] =\n",
       " CartesianIndex(1, 2, 4)  CartesianIndex(1, 3, 4)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I1 = CartesianIndex(1, 2, 2)\n",
    "I2 = CartesianIndex(1, 3, 4)\n",
    "\n",
    "I1:I2 # CartesianIndices 实际上就是 Array{CartesianIndex}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×2×3 Array{Int64,3}:\n",
       "[:, :, 1] =\n",
       " 9  11\n",
       "\n",
       "[:, :, 2] =\n",
       " 15  17\n",
       "\n",
       "[:, :, 3] =\n",
       " 21  23"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[I1:I2] # 等价于 A[1:1, 2:3, 2:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "也可以用 `CartesianIndices` 来快速构造对应矩阵的下标"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3×4 CartesianIndices{3,Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}}:\n",
       "[:, :, 1] =\n",
       " CartesianIndex(1, 1, 1)  CartesianIndex(1, 2, 1)  CartesianIndex(1, 3, 1)\n",
       " CartesianIndex(2, 1, 1)  CartesianIndex(2, 2, 1)  CartesianIndex(2, 3, 1)\n",
       "\n",
       "[:, :, 2] =\n",
       " CartesianIndex(1, 1, 2)  CartesianIndex(1, 2, 2)  CartesianIndex(1, 3, 2)\n",
       " CartesianIndex(2, 1, 2)  CartesianIndex(2, 2, 2)  CartesianIndex(2, 3, 2)\n",
       "\n",
       "[:, :, 3] =\n",
       " CartesianIndex(1, 1, 3)  CartesianIndex(1, 2, 3)  CartesianIndex(1, 3, 3)\n",
       " CartesianIndex(2, 1, 3)  CartesianIndex(2, 2, 3)  CartesianIndex(2, 3, 3)\n",
       "\n",
       "[:, :, 4] =\n",
       " CartesianIndex(1, 1, 4)  CartesianIndex(1, 2, 4)  CartesianIndex(1, 3, 4)\n",
       " CartesianIndex(2, 1, 4)  CartesianIndex(2, 2, 4)  CartesianIndex(2, 3, 4)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CartesianIndices(A) # A 的所有下标"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以及坐标的偏移等计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:\n",
       " CartesianIndex(2, 2)  CartesianIndex(2, 3)  CartesianIndex(2, 4)\n",
       " CartesianIndex(3, 2)  CartesianIndex(3, 3)  CartesianIndex(3, 4)\n",
       " CartesianIndex(4, 2)  CartesianIndex(4, 3)  CartesianIndex(4, 4)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "δ = CartesianIndex(1, 1)\n",
    "I = CartesianIndex(3, 3)\n",
    "\n",
    "# 构造了一组以(3, 3)为中心的坐标\n",
    "I-δ:I+δ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:\n",
       " CartesianIndex(0, 0)  CartesianIndex(0, 1)  CartesianIndex(0, 2)\n",
       " CartesianIndex(1, 0)  CartesianIndex(1, 1)  CartesianIndex(1, 2)\n",
       " CartesianIndex(2, 0)  CartesianIndex(2, 1)  CartesianIndex(2, 2)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I_first = CartesianIndex(1, 1)\n",
    "I_last = CartesianIndex(5, 5)\n",
    "\n",
    "I = CartesianIndex(1, 1)\n",
    "δ = CartesianIndex(1, 1)\n",
    "\n",
    "I-δ:I+δ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了让 `I-δ:I+δ` 在有效的区域内，比如说在 `I_first:I_last` 之间，可以使用 `min`, `max` 等工具来进行处理\n",
    "\n",
    "以下两种都是常见的处理策略："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 CartesianIndices{2,Tuple{UnitRange{Int64},UnitRange{Int64}}}:\n",
       " CartesianIndex(1, 1)  CartesianIndex(1, 2)\n",
       " CartesianIndex(2, 1)  CartesianIndex(2, 2)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(I_first, I-δ) : min(I_last, I+δ)\n",
    "\n",
    "# 一些负的坐标被移除了"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "大致理解的话，`CartesianIndex` 可以帮助我们将多维下标的处理抽象出来，从而来避免复杂的for循环\n",
    "\n",
    "回到之前的均值滤波操作，我们可以用 `CartesianIndex` 将它改写成一般形式\n",
    "\n",
    "```julia\n",
    "function mean_filter_1d(X::AbstractVector; window_size=3)\n",
    "    out = similar(X)\n",
    "    \n",
    "    δ = window_size ÷ 2\n",
    "    for i in 1:length(X)\n",
    "        # min, max确保下标不溢出\n",
    "        l = max(1, i-δ)\n",
    "        r = min(length(X), i+δ)\n",
    "        \n",
    "        # 计算一个局部的均值\n",
    "        out[i] = mean(X[l:r])\n",
    "    end\n",
    "    return out\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mean_filter (generic function with 1 method)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 这个函数可以对任意维的矩阵进行迭代\n",
    "\n",
    "function mean_filter(X; window_size=3)\n",
    "    out = similar(X)\n",
    "    \n",
    "    δ = CartesianIndex(ntuple(_ -> window_size÷2, ndims(X)))\n",
    "\n",
    "    R = CartesianIndices(X)\n",
    "    I_first, I_last = first(R), last(R)\n",
    "    \n",
    "    for i in R\n",
    "        Rᵢ = max(I_first, i-δ) : min(I_last, i+δ)\n",
    "        out[i] = mean(X[Rᵢ])\n",
    "    end\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 一维\n",
    "X = Float64.(rand(1:5, 10))\n",
    "mean_filter(X);\n",
    "\n",
    "# 二维\n",
    "X = Float64.(rand(1:5, 10, 10))\n",
    "mean_filter(X);\n",
    "\n",
    "# 3维 -- 滑动窗口是一个三维立体\n",
    "X = Float64.(rand(1:5, 10, 10, 3))\n",
    "mean_filter(X);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. view and copy\n",
    "\n",
    "在对矩阵进行取下标操作时，默认的方式是 `copy`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3×4 Array{Int64,3}:\n",
       "[:, :, 1] =\n",
       " 1  3  5\n",
       " 2  4  6\n",
       "\n",
       "[:, :, 2] =\n",
       " 7   9  11\n",
       " 8  10  12\n",
       "\n",
       "[:, :, 3] =\n",
       " 13  15  17\n",
       " 14  16  18\n",
       "\n",
       "[:, :, 4] =\n",
       " 19  21  23\n",
       " 20  22  24"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Int64,2}:\n",
       " 0  0  0\n",
       " 0  0  0"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s1 = A[:, :, 1]\n",
    "\n",
    "fill!(s1, 0) # 因为 s1 是通过复制得到的，所以修改 s1 不会修改 A 的值\n",
    "s1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3×4 Array{Int64,3}:\n",
       "[:, :, 1] =\n",
       " 1  3  5\n",
       " 2  4  6\n",
       "\n",
       "[:, :, 2] =\n",
       " 7   9  11\n",
       " 8  10  12\n",
       "\n",
       "[:, :, 3] =\n",
       " 13  15  17\n",
       " 14  16  18\n",
       "\n",
       "[:, :, 4] =\n",
       " 19  21  23\n",
       " 20  22  24"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这样做是比较安全的，可以避免以外的数据修改，但在很多时候，这样做会带来大量额外的内存开销。\n",
    "\n",
    "在 Julia 下，可以通过 `view`, `@view` 或者 `@views` 来避免数据的复制"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 view(::Array{Int64,3}, :, :, 1) with eltype Int64:\n",
       " 0  0  0\n",
       " 0  0  0"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s2 = @view A[:, :, 1] # view(A, :, :, 1)\n",
    "fill!(s2, 0) # 因为 s1 是通过 view 得到的，所以修改 s1 会同样修改 A 的值\n",
    "s2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3×4 Array{Int64,3}:\n",
       "[:, :, 1] =\n",
       " 0  0  0\n",
       " 0  0  0\n",
       "\n",
       "[:, :, 2] =\n",
       " 7   9  11\n",
       " 8  10  12\n",
       "\n",
       "[:, :, 3] =\n",
       " 13  15  17\n",
       " 14  16  18\n",
       "\n",
       "[:, :, 4] =\n",
       " 19  21  23\n",
       " 20  22  24"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之前的 `mean_filter` 这个函数就可以通过这一点来进行优化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  24.002 μs (303 allocations: 69.38 KiB)\n"
     ]
    }
   ],
   "source": [
    "# 这个函数可以对任意维的矩阵进行迭代\n",
    "\n",
    "function mean_filter(X; window_size=3)\n",
    "    out = similar(X)\n",
    "    \n",
    "    δ = CartesianIndex(ntuple(_ -> window_size÷2, ndims(X)))\n",
    "\n",
    "    R = CartesianIndices(X)\n",
    "    I_first, I_last = first(R), last(R)\n",
    "    \n",
    "    for i in R\n",
    "        Rᵢ = max(I_first, i-δ) : min(I_last, i+δ)\n",
    "        out[i] = mean(X[Rᵢ])\n",
    "    end\n",
    "    return out\n",
    "end\n",
    "\n",
    "@btime mean_filter($X);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  14.130 μs (3 allocations: 2.56 KiB)\n"
     ]
    }
   ],
   "source": [
    "# 这个函数可以对任意维的矩阵进行迭代\n",
    "\n",
    "function mean_filter(X; window_size=3)\n",
    "    out = similar(X)\n",
    "    \n",
    "    δ = CartesianIndex(ntuple(_ -> window_size÷2, ndims(X)))\n",
    "\n",
    "    R = CartesianIndices(X)\n",
    "    I_first, I_last = first(R), last(R)\n",
    "    \n",
    "    for i in R\n",
    "        Rᵢ = max(I_first, i-δ) : min(I_last, i+δ)\n",
    "        out[i] = mean(@view X[Rᵢ]) # 这里只创建 view 从而避免额外的性能开销\n",
    "    end\n",
    "    return out\n",
    "end\n",
    "\n",
    "@btime mean_filter($X);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.0",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
